clear; clc; close all;

% 迭代总步数及舍去的前期步数(瞬态)
total_iter = 10000;  
transient  = 500;   
record_num = total_iter - transient;  % 实际记录点数


%% (a) 2D-MLM
mu_mlm = 0.2; 
k_mlm  = 2.08;
x0_mlm = 0.2; 
q0_mlm = 0.1;
[xMLM, qMLM] = attractor(@MLM, mu_mlm, k_mlm, x0_mlm, q0_mlm, total_iter, transient);


%% (b) 2D-MSM
mu_msm = -0.1; 
k_msm  = 1.88;
x0_msm = 0.8; 
q0_msm = 0.6;
[xMSM, qMSM] = attractor(@MSM, mu_msm, k_msm, x0_msm, q0_msm, total_iter, transient);


%% (c) 2D-MTM
mu_mtm = -0.6; 
k_mtm  = 1.78;
x0_mtm = 0.8; 
q0_mtm = -0.8;
[xMTM, qMTM] = attractor(@MTM, mu_mtm, k_mtm, x0_mtm, q0_mtm, total_iter, transient);


%% (d) 2D-SMM
mu_smm = 0.5; 
k_smm  = 2.33;
x0_smm = 1; 
q0_smm = 1;
[xSMM, qSMM] = attractor(@SMM, mu_smm, k_smm, x0_smm, q0_smm, total_iter, transient);


%% 绘图
% (a) 2D-MLM
subplot(2,2,1);
plot(xMLM, qMLM, 'm.', 'MarkerSize', 2);
xlabel('x'); ylabel('q');
title('2D-MLM (\mu=0.2, k=2.08)');
axis equal; grid on;  % 让坐标比例一致，便于看相图形状

% (b) 2D-MSM
subplot(2,2,2);
plot(xMSM, qMSM, 'g.', 'MarkerSize', 2);
xlabel('x'); ylabel('q');
title('2D-MSM (\mu=-0.1, k=1.88)');
axis equal; grid on;

% (c) 2D-MTM
subplot(2,2,3);
plot(xMTM, qMTM, 'r.', 'MarkerSize', 2);
xlabel('x'); ylabel('q');
title('2D-MTM (\mu=-0.6, k=1.78)');
axis equal; grid on;

% (d) 2D-SMM
subplot(2,2,4);
plot(xSMM, qSMM, 'k.', 'MarkerSize', 2);
xlabel('x'); ylabel('q');
title('2D-SMM (\mu=0.5, k=2.33)');
axis equal; grid on;


%% 生成吸引子
function [xdata, qdata] = attractor(map_func, mu, k, x0, q0, total_iter, transient)
    x = x0;
    q = q0;
    
    % 先做 transient 步以排除瞬态
    for i = 1:transient
        [x, q] = map_func(x, q, mu, k);
    end
    
    % 再记录后续轨迹
    record_num = total_iter - transient;
    xdata = zeros(1, record_num);
    qdata = zeros(1, record_num);
    
    xdata(1) = x;
    qdata(1) = q;
    for i = 2:record_num
        [xdata(i), qdata(i)] = map_func(xdata(i-1), qdata(i-1), mu, k);
    end
end


%% 1) 2D-MLM 映射
function [x_next, q_next] = MLM(x, q, mu, k)
    % x_{n+1} = mu*x_n*(1 - x_n) + k*cos(q_n)*x_n
    % q_{n+1} = q_n + x_n
    x_next = mu*x*(1 - x) + k*cos(q)*x;
    q_next = q + x;
end

%% 2) 2D-MSM 映射
function [x_next, q_next] = MSM(x, q, mu, k)
    % x_{n+1} = mu*sin(2*pi*x_n) + k*cos(q_n)*x_n
    % q_{n+1} = q_n + x_n
    x_next = mu*sin(2*pi*x) + k*cos(q)*x;
    q_next = q + x;
end

%% 3) 2D-MTM 映射
function [x_next, q_next] = MTM(x, q, mu, k)
    % 分段
    if x < 0.5
        x_next = mu*x + k*cos(q)*x;
    else
        x_next = mu*(1 - x) + k*cos(q)*x;
    end
    q_next = q + x;
end

%% 4) 2D-SMM 映射
function [x_next, q_next] = SMM(x, q, mu, k)
    % x_{n+1} = mu*x_n + k*cos(q_n)*x_n
    % q_{n+1} = q_n + x_n
    x_next = mu*x + k*cos(q)*x;
    q_next = q + x;
end
